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Abstract. LISA (Laser Interferometer Space Antenna) is a proposed space 
mission, which will use coherent laser beams exchanged between three remote 
spacecraft to detect and study low-frequency cosmic gravitational radiation. In 
the low-part of its frequency band, the LISA strain sensitivity will be dominated 
by the incoherent superposition of hundreds of millions of gravitational wave 
signals radiated by inspiraling white-dwarf binaries present in our own galaxy. 
In order to estimate the magnitude of the LISA response to this background, we 
have simulated a synthesized population that recently appeared in the literature. 
Our approach relies on entirely analytic expressions of the LISA Time-Delay 
Interferometric responses to the gravitational radiation emitted by such systems, 
which allows us to implement a computationally efficient and accurate simulation 
of the background in the LISA data. We find the amplitude of the galactic white- 
dwarf binary background in the LISA data to be modulated in time, reaching 
a minimum equal to about twice that of the LISA noise for a period of about 
two months around the time when the Sun-LISA direction is roughly oriented 
towards the Autumn equinox. This suggests that, during this time period, LISA 
could search for other gravitational wave signals incoming from directions that 
are away from the galactic plane. Since the galactic white-dwarfs background will 
be observed by LISA not as a stationary but rather as a cyclostationary random 
process with a period of one year, we summarize the theory of cyclostationary 
random processes and present the corresponding generalized spectral method 
needed to characterize such process. We find that, by measuring the generalized 
spectral components of the white-dwarf background, LISA will be able to infer 
properties of the distribution of the white-dwarfs binary systems present in our 
Galaxy. 
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1. Introduction 

The Laser Interferometric Space Antenna (LISA) is a space mission jointly proposed 
to the National Aeronautics and Space Administration (NASA) and the European 
Space Agency (ESA). Its aim is to detect and study gravitational waves (GW) in 
the millihertz frequency band. It will use coherent laser beams exchanged between 
three identical spacecraft forming a giant (almost) equilateral triangle of side 5 x 10 6 
kilometers. By monitoring the relative phase changes of the light beams exchanged 
between the spacecraft, it will extract the information about the gravitational waves 
it will observe at unprecedented sensitivities pQ. 

The astrophysical sources that LISA is expected to observe within its operational 
frequency band (10~ 4 — 1 Hz) include extra-galactic super-massive black-hole 
coalescing binaries, stochastic gravitational wave background from the early universe, 
and galactic and extra-galactic coalescing binary systems containing white dwarfs and 
neutron stars. 

Recent surveys have uniquely identified twenty binary systems emitting 
gravitational radiation within the LISA band, while population studies have concluded 
that the large number of binaries present in our own galaxy should produce a stochastic 
background that will lie significantly above the LISA instrumental noise in the low-part 
of its frequency band. It has been shown in the literature (see for a recent study 
and 0] for earlier investigations) that these sources will be dominated by detached 
white-dwarf — white-dwarf (WD- WD) binaries, with 1.1 x 10 s of such systems in our 
Galaxy. By using the distribution for the parameters of the WD- WD binaries given 
in (J3j) we have simulated the LISA response to the WD-WD background. 

This paper is organized as follows. In Section[21we recall the analytic expression of 
the unequal-arm Michelson combination, X, of the LISA Time-Delay Interferometric 
(TDI) response to a signal radiated by a binary system. In Section [3| we give a 
summary of how the WD-WD binary population was obtained, and a description 
of our numerical simulation of the X response to it. In Section 01 we describe the 
numerical implementation of our simulation of the LISA X response to the WD-WD 
background, and summarize our results. The time-dependence and periodicity of the 
magnitude of the WD-WD galactic background in the LISA data implies that it is 
not a stationary but rather a cyclostationary random process of period one year. In 
Section[S]we provide a brief summary of the theory of cyclostationary random processes 
relevant to the LISA detection of the WD-WD galactic background, and apply it to 
three years worth of simulated LISA X data in Sectional We find that, by measuring 
the generalized spectral components of such cyclostationary random process, LISA 
will be able to infer key-properties of the distribution of the WD-WD binary systems 
present in our own Galaxy. 

The non-stationarity of the WD-WD background was first pointed out 
by Giampieri and Polnarev j!8l under the assumption of sources distributed 
anisotropically, and they also obtained the Fourier expansion of the sample variance 
and calculated the Fourier coefficient for simplified WD-WD binary distributions in the 
Galactic disc. What was however not realized in their work is that this non-stationary 
random process is actually cyclostationary, i.e. there exists cyclic spectra that can in 
principle allow us to infer more information about the WD-WD background than one 
could obtain by just estimating the zero-order spectrum. 
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2. The LISA response to signals from binary systems 



The gravitational wave response X GW (t) of the first generation Michelson TDI 
combination to a signal from a binary system has been derived in and it can 
be written in the following form 



X GW (t) = Re \A(x,t) e-^ (t) 



(1) 



where x — uj s L (lu s being the angular frequency of the GW signal in the source 
reference frame), and the expressions for the complex amplitude A(x,t) and the real 
phase <j)(t) are 



A(x, t) — 2 x sm(ir) j 
— sinc[(l — 



sinc[(l + ca(t))|] e ix (i+ d2 «) + smc[(l - c 2 (i))|] 
<*(*))f] ^ 



a ixC±+d 2 (t)) 



<j){t) = us s t + uj s R cos (3 cos(fii + ?/o — A) . 

In equation J3J i? is the distance of the guiding center of the LISA array, o, from 
the Solar System Barycenter (SSB), (/3, A) are the ecliptic latitude and longitude 
respectively of the source location in the sky, = 27r/year, and ijq defines the position 
of the LISA guiding center in the ecliptic plane at time t — 0. Note that the functions 
Ck(t), dk{t), and Bk{t) (k — 2,3) do not depend on x. The analytic expressions for 
Ck{t), and dk(t) are the same as those given in equations (46,47) of reference |§J, while 
the functions Bk{t) (k = 2,3) are equal to 



B k (t) = (a« + i a < 3 ») Ufc (t) + (a (2) + i a (4) ) v k (t) 



(4) 



The coefficients (a^',a 



(2) a (3) „.(4) 



a 1 



depend only on the two independent amplitudes 
of the gravitational wave signal, (h + , h x ), the polarization angle, ip, and an arbitrary 
phase, that the signal has at time t = 0. Their analytic expressions are given in 
equations (41-44) of reference [S], while the functions Uk(t), and Vk(t) (k = 2,3) are 
given in equations (27,28) in the same reference. 

Since most of the gravitational wave energy radiated by the galactic WD- WD 
binaries will be present in the lower part of the LISA sensitivity frequency band, say 
between 10~ 4 — 10~ 3 Hz, it is useful to provide an expression for the Taylor expansion 
of the X response in the long- wavelength limit (LWL), i.e. when the wavelength 
of the gravitational wave signal is much larger than the LISA armlength (x « 1). 
The LWL expression will allow us to analytically describe the general features of the 
white dwarfs background in the A-combination, and derive computationally efficient 
algorithms for numerically simulating the WD- WD background in the LISA data. 



The nth-order truncation, X 1 ^ (t), of the Taylor expansion of X GW 
series of x can be written in the following form 



X^(t) = ReJ2A( k \t) 



x k+2 e 



-i<p(t) 



(t) in power 



(5) 



fe=0 



where the first three functions of time A^ k \t), k < 2 are equal to 
A<°> = 4 [B 2 - B 3 ] , 
A« = 4z [(dj, + 2) B 2 - + 2)B 3 ] , 
28 1 
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Note that the form we adopted for X GW (t) (equation ^) makes the derivation of the 
functions (t) particularly easy since the dependence on x in A(x, t) is now limited 
only to the coefficients in front of the two functions B2 (i) and B3 (t) (see equation @ ) . 

Although it is generally believed that the lowest order long-wavelength expansion 
of the X combination, X ^ , is sufficiently accurate in representing a gravitational 
wave signal in the low-part of the LISA frequency band, there has not been in the 
literature any quantitative analysis of the error introduced by relying on such a zero- 
order approximation. Since any TDI combination will contain a linear superposition 
of tens of millions of signals, it is crucial to estimate such an error as a function of 
the order of the approximation, n. In order to determine how many terms we needed 
to use for a given signal angular frequency, uj s , we relied on the following " matching 
function" 



Af(A GW ,A ( ^- 



J^[X^(t)-X G J(t)] 2 dt 
\ tf{X™(t)] 2 dt 

Equation J7J estimates the percent root-mean-squared error implied by using the n th 
order LWL approximation. With T — lyear we have found [7] that at 5 x 10~ 4 
Hz, for instance, the zero-order LWL approximation (n = 0) of the X combination 
shows an r.m.s. deviation from the exact response equal to about 10 percent. As 
expected, this inaccuracy increases for signals of higher frequencies, becoming equal 
to 40 percent at 2 x 10 -3 Hz. With n = 1 the accuracy improves showing that 
the X£W response deviates from the exact one with an r.m.s. error smaller than 10 
percent in the frequency band (10~ 4 — 2 x 10~ 3 ) Hz. In our simulation we have actually 
implemented the n = 2 LWL expansion because it was possible and easy to do. 



3. White Dwarf binary population distribution 

The gravitational wave signal radiated by a WD- WD binary system depends on 
eight parameters, (0 O , l, ip, D, /?, A, M. c , w s ), which are the constant phase of the signal 
(4> ) at the starting time of the observation, the inclination angle (t) of the angular 
momentum of the binary system relative to the line of sight, the polarization angle (ip) 
describing the orientation of the wave polarization axes, the distance (D) to the binary, 
the angles (A, (3) describing the location of the source in the sky relative to the ecliptic 
plane, the chirp mass (Mc), and the angular frequency (uj s ) in the source reference 
frame respectively. Since it can safely be assumed that the chirp mass M c and the 
angular frequency uj s are independent of the source location and of the remaining 
angular parameters <p Q , l, tp, and because there are no physical arguments for preferred 
values of the constant phase <p and the orientation of the binary given by the angles 
t and ip, it follows that the joint probability distribution, P((p , t, tp, D, (3, A, M c , Us), 
can be rewritten in the following form 

P{ct> , l, i>, D, /3, A, Mc, w.) = P x {<t>o)P2 (O-Pa (V0 p 4 (D, (3, X)P 5 (M C , u„) . 
In the implementation of our simulation we have assumed the angles <j) and tp to 
be uniformly distributed in the interval [0, 2-7r), and cost uniformly distributed in the 
interval [—1,1]. We further assumed the binary systems to be randomly distributed 
in the Galactic disc according to the following axially symmetric distribution Va{R, z) 
(see [12] Eq. (5)) 

V ^ z >~ ^Jf2 - (9) 
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where (-R, z) are cylindrical coordinates with origin at the galactic center, H — 2.5 
kpc, and z a — 200 pc, and it is proportional to Pa{D, A, [3) through the Jacobian of the 
coordinate transformation. Note that the position of the Sun in this coordinate system 
is given by Rq = 8.5 kpc and Zq = —30 pc. We then generate the positions of the 
sources from the distribution given by Eq. @ and map them to their corresponding 
ecliptic coordinates (D,f3,X). 

The physical properties of the WD- WD population (M c = (mim2) 3 ' 5 / (mi + TO2) 1 
with mi, 7712 being the masses of the two stars, and lo s — 2?:f s = 47r/P or b) are taken 
from the binary population synthesis simulation discussed in |14j . For details on this 
simulation we r efer the reader to [H], and for earlier work to |3l II ITHl ITT1 IT2l IH] . 
The basic ingredient for these simulations is an approximate binary evolution code. A 
representation of the complete Galactic population of binaries is produced by evolving 
a large (typically 10 6 ) number of binaries from their formation to the current time, 
where the distributions of the masses and separations of the initial binaries are esti- 
mated from the observed properties of local binaries. This initial-to-final parameter 
mapping is then convolved with an estimate of the binary formation rate in the history 
of the Galaxy to obtain the total Galactic population of binaries at the present time. 
From these the binaries of interest can then be selected. In principle this technique is 
very powerful, although the results can be limited by the limited knowledge we have 
on many aspects of binary evolution. For WD- WD binaries, the situation is better 
than for many other populations, since the observed population of WD- WD binaries 
allows us to gauge the models (e.g. [T^]'). 

We also include the population of semi-detached WD-WD binaries (usually 
referred to as AM CVn systems) that are discussed in detail in [2}. In these binaries 
one white dwarf transfers its outer layers onto a companion white dwarf. Due to the 
redistribution of mass in the system, the orbital period of these binaries increases in 
time, even though the angular momentum of the binary orbit still decreases due to 
gravitational wave losses. The formation of these systems is very uncertain, mainly 
due to questions concerning the stability of the mass transfer (e.g. [TB)1. 

From the models of the Galactic population of the detached WD-WD binaries 
and AM CVn systems two dimensional histograms were created, giving the expected 
number of both WD-WD binaries and AM CVn systems currently present in the 
Galaxy as function of the log of the GW radiation frequency, f s {= ^ s /27r) and chirp 
mass, M. c . In the case of the detached WD-WD binaries, the (log/ s , Ai c ) space was 
defined over the set log/ s € [—6,-1], M. c E (0,1.5] and contained 50 x 30 grid 
points, while in the case of the AM CVn systems the region is intrinsically smaller, 
log/, e [—4,-1.5], M. c £ (0,1.2], containing only 25 x 24 grid points. In order 
to generate values for these distributions within each grid rectangle we have used an 
interpolation method 0. 

Figure JQ) shows the distribution of the number of detached WD-WD binaries as 
a function of the signal frequency and chirp mass in the form of a contour plot. This 
distribution reaches its maximum within the LISA frequency band when the chirp 
mass is equal to ~ 0.25 A4q, and it monotonically decreases as a function of the 
signal frequency. The distribution of the number of AM CVn systems has instead a 
rather different shape, as shown by the contour plot given in Figure 
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Figure 1. Distribution of detached white-dwarf — white-dwarf binaries in our 
galaxy as a function of the gravitational wave frequency and chirp mass 
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Figure 2. Distribution of AM CVn binary systems in our galaxy as a function 
of the gravitational wave frequency and chirp mass 
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4. Numerical Simulation 

In order to simulate the LISA X response to the population of WD- WD binaries 
derived in Section|21one needs to coherently add the LISA response to each individual 
signal. Although this could naturally be done in the time domain, the actual CPU 
time required to successfully perform such a simulation would be unacceptably long. 
The generation in the time domain of one year of X GW (t) response to a single signal 
sampled at a rate of 16 seconds would require about 1 second with an optimized CH — h 
code running on a Pentium IV 3.2 GHz processor. Since the number of signals from 
the background is of the order 10 8 , it is clear that a different algorithm is needed 
for simulating the background in the LISA data within a reasonable amount of time. 
We were able to derive and implement numerically an analytic formula of the Fourier 
transform of each binary signal, which has allowed us to reduce the computational 
time by almost a factor 100. 

We have obtained an approximate analytic expression for the finite-time Fourier 
transform of each WD- WD signal using the LWL expansion (JSJ with n = 2 and 
by applying the Nuttall's modified Blackman-Harris window 19- in order to avoid 
spectral leakage. We have then coherently added in the Fourier domain all the signals 
radiated by the WD- WD galactic binary population described in section ©. After 
inverse Fourier transforming the synthesized response and then removing the window 
from it, we finally obtained the time-domain representation of the background as it 
will be seen in the LISA TDI combination X. This is shown in Figure ©, where we 
plot three years worth of simulated X GW (t), and include the LISA noise pQ. The one- 
year periodicity induced by the motion of LISA around the Sun is clearly noticeable. 
One other interesting feature shown by Figure © is that the amplitude response 
reaches absolute minima when the Sun-LISA direction is roughly oriented towards the 
Autumn equinox, while the absolute maxima take place when the Sun-LISA direction 
is oriented roughly towards the Galactic center This is because the ecliptic plane 
is not parallel to the galactic plane, and our own solar system is about 8.5 kpc away 
from the galactic center (where most of the of WD- WD binaries are concentrated). 
As a result it follows that the LISA X GW response to the WD- WD background does 
not have a six-months periodicity. 

Note also that, for a time period of about 2 months, the absolute minima reached 
by the amplitude of the LISA response to the WD- WD background is only a factor less 
than 2 larger than the level of the instrumental noise. This implies that during these 
observation times LISA should be able to search for other sources of gravitational 
radiation that are not located in the galactic plane. This might turn out to be 
the easiest way to mitigate the detrimental effects of the WD- WD background when 
searching for other sources of gravitational radiation. We will quantitatively analyze 
in a follow up work how to take advantage of this observation in order to optimally 
search, during these time periods, for sources that are off the galactic plane. 

5. Cyclostationary processes 

The results of our simulation (Figure|3| indicate that the LISA X GW response to the 
background can be regarded, in a statistical sense, as a periodic function of time. This 
is consequence of the deterministic (and periodic) motion of the LISA array around the 
Sun. Since its autocorrelation function will also be a periodic function of period one 
year, it follows that any LISA response to the WD- WD background should no longer be 
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Figure 3. Simulated time series of the WD-WD Galactic background signal 
of 3 years of data. The time series of LISA instrumental noise is displayed for 
comparison. 

treated as a stationary random process but rather as a periodically correlated random 
process. These kind of processes have been studied for many years, and are usually 
referred to as cyclostationary random processes (see |2(J) for a comprehensive overview 
of the subject and for more references). In what follows we will briefly summarize the 
properties of cyclostationary processes that are relevant to our problem. 

A continuous stochastic process X(t) having finite second order moments is said 
to be cyclostationary with period T if the following expectation values 

E[X{t)\ =m(t) =m(t + T), (10) 

E[X{t')X(t)] = C(t', t) = C{t' + T,t + T) (11) 

are periodic functions of period T, for every (t',t) € R x R. For simplicity from now 
on we will assume m(t) = 0. 

If X(t) is cyclostationary, then the function B(t, r) = C(t + r, t) for a given r S R 
is periodic with period T, and it can be represented by the following Fourier series 

oo 

B(t,r)= £ B r (ry 2 ^ , (12) 
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where the functions B r (r) are given by 



T 



B r (r) = ~ J B(t,r)e-^ r * dt . (13) 

The Fourier transforms g r (f) of B r (r) are the so called "cyclic spectra" of the 
cyclostationary process X(t) [2U| 

/oo 
B r (T)e~^ dr . (14) 
-oo 

If a cyclostationary process is real, the following relationships between the cyclic 
spectra hold 

B_ r (r) =B;(t) , (15) 
9-A-f) = 9* r (f) , (16) 
where the symbol * means complex conjugation. This implies that, for a real 
cyclostationary process, the cyclic spectra with r > contain all the information 
needed to characterize the process itself. 

The function <t 2 (t) = 5(0, r) is the variance of the cyclostationary process X(t), 
and it can be written as a Fourier decomposition as a consequence of Eq. (|13fl 

oo 

o\t)= H r ^, (17) 

r— — oo 

where H r = B r (0) are harmonics of the variance a 2 . From Eq. 11511 it follows that 
H_ r = H*. 

For a discrete, finite, real time series X t , t = 1, . . . , N we can estimate the cyclic 
spectra by generalizing standard methods of spectrum estimation used with stationary 
processes. Assuming again the mean value of the time series X t to be zero, the cyclic 
autocorrelation sequences are defined as 

1 N ~ W 

s i = nJ2 x t x l+ \i\ e ~ !Ssr ^ ■ ( 18 ) 

t=l 

It has been shown [20] that the cyclic autocorrelations are asymptotically (i.e. for 
N — > oo) unbiased estimators of the functions B r (r). The Fourier transforms of the 
cyclic autocorrelation sequences s\ are estimators of the cyclic spectra g r {f)- These 
estimators are asymptotically unbiased, and are called "inconsistent estimators" of 
the cyclic spectra, i.e. their variances do not tend to zero asymptotically. In the case 
of Gaussian processes [201 consistent estimators can be obtained by first applying a 
lag window to the cyclic autocorrelation and then perform a Fourier transform. This 
procedure represents a generalization of the well-known technique for estimating the 
spectra of stationary random processes |21| . 

An alternative procedure for identifying consistent estimators of the cyclic spectra 
is to first take the Fourier transform, X(f), of the time series X(t) 

N 

x(f) =Y, x t e ~ i2nnt ' 1) ( 19 ) 

and then estimate the cyclic periodograms g r (f) 

X(f)X*(f-^) 
9r(f) = Jj ■ (20) 
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By finally smoothing the cyclic periodograms, consistent estimators of the spectra 
g r (f) are then obtained. The estimators of the harmonics H r of the variance a 2 of 
a cyclostationary random process can be obtained by first forming a sample variance 
of the time series X t . The sample variance is obtained by dividing the time series 
X t into contiguous segments of length tq such that To is much smaller than the 
period T of the cyclostationary process, and by calculating the variance erf over 
each segment. Estimators of the harmonics can be obtained by Fourier analyzing 
the series erf. Note that the definitions of (i) zero order (r = 0) cyclic autocorrelation, 
(ii) periodogram, and (iii) zero order harmonic of the variance, coincide with those 
usually adopted for stationary random processes. Thus, even though a cyclostationary 
time series is not stationary, ordinary spectral analysis can be used for obtaining the 
zero order spectrum. Note, however, that cyclostationary random processes provide 
more spectral information about the time series they are associated with due to the 
existence of cyclic spectra with r > 0. 

As an important and practical application, let us consider a time series yt 
consisting of the sum of a stationary random process, n t , and a cyclostationary one X t 
(i.e. yt = n t +X t ). Let the variance of the stationary time series be v 1 and its spectral 
density be £ (/). It is easy to see that the resulting process is also cyclostationary. If 
the two processes are uncorrelatcd, then the zero order harmonic Sq of the variance 
of the combined processes is equal to 

E 2 = u 2 + a 2 , (21) 

and the zero order spectrum Go(/) of yt is 

G (f) =£(/)+; 9o (/) • (22) 

The harmonics of the variance as well as the cyclic spectra of yt with r > coincide 
instead with those of X t . In other words, the harmonics of the variance and the cyclic 
spectra of the process yt with r > contain information only about the cyclostationary 
process X t , and are not "contaminated" by the stationary process nt- 

6. Data analysis of the background signal 

We have numerically implemented the methods outlined in Section [S] and applied 
them to our simulated WD- WD background signal. A comparison of the results 
of our simulation of the detached WD- WD background with the calculation of the 
background by Hils and Bender ^3 El is shown in Figure^] Note that in Figure0]we 
have given the spectra of the background that do not include LISA instrumental noise. 
We find that the amplitude of the background from our simulation is a factor of more 
than 2 smaller than that of Hils and Bender. The level of the WD- WD background 
is determined by the number of such systems in the Galaxy. We estimate that our 
number WD- WD binaries should be correct within a factor 5 and thus the amplitude of 
the background should be right within a factor of In Figure0]we have plotted the 
two backgrounds against the LISA spectral density and we have also included the LISA 
sensitivity curve. The latter is obtained by dividing the instrumental noise spectral 
density by the detector GW transfer function averaged over isotropically distributed 
and randomly polarized signals. In the zero-order long wavelength approximation this 
averaged transfer function is equal to y3/20. Note that for our simulation in the 
region of the LISA band below 0.2 millihertz the power of the WD- WD background 
is smaller than that of the instrumental noise. 
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Figure 4. Comparison of detached WD-WD background obtained from 
binary population synthesis simulation ( I '■'> 1141 1 with the WD-WD background 
calculated by Hils and Bender Amplitude spectral density of LISA 

instrumental noise and LISA sensitivity curve are drawn for comparison. All 
spectral densities are one-sided. 



Our analysis was applied to 3 years of LISA X data consisting of a coherent 
superposition of signals emitted by detached WD-WD binaries, by semi-detached 
binaries (AM CVn systems), and of simulated instrumental noise. The noise was 
numerically generated by using the spectral density of the TDI X observable given in 
HJ|. In addition a 1 mHz low-pass filter was applied to our data set in order to focus 
our analysis to the frequency region in which the WD-WD stochastic background is 
expected to be dominant. 

The results of the Fourier analysis of the sample variance of the background 
signal are shown in Figure [3] The top panel of Figure [S] shows the sample variance 
of the simulated data for which the variances were estimated over a period of 1 week; 
periodicity is clearly visible. The bottom panel instead shows the Fourier analysis of 
the sample variance for which we have removed the mean from the sample variance 
time series. The vertical lines correspond to multiples of 1 year; two harmonics can 
clearly be distinguished from noise. The other peaks of the spectrum that fall roughly 
half way between the multiples of 1 /year frequency, are from the rectangular window 
inherent to the finite time series. In our analysis we have plotted the absolute values 
\H r \ of the coefficients of the Fourier decomposition given by En.l|17|l. By using the 
Fourier transform method, we cannot resolve higher than the second harmonic in our 
3-year data set. We conclude that only the two first harmonics can be extracted 
reliably from the data. 

It is useful to compare the results of our numerical analysis against the analytic 
calculations of Giampieri and Polnarev ^Hj- Their analytic expressions for the 
harmonics of the variance of a background due to binary systems distributed in the 



Simulation of the White Dwarf — White Dwarf galactic background in the LISA data.Yl 




Fourier analysis of the variance 
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Figure 5. Top panel: sample variance of the simulated WD-WD background. 
Three years of data are simulated. Data include two populations of the WD-WD 
binaries, detached and semi-detached ones added to the LISA instrumental noise. 
The data is passed through a low-pass filter with a cut-off frequency of 1 mHz. 
Bottom panel: Fourier analysis of the sample variance. Two harmonic are clearly 
resolved. 

galactic disc are given in Eq. (42) and shown in Figure 4 of Our estimation 

roughly matches theirs in that the 0th order harmonics is dominant and the first two 
harmonics have more power than the remaining ones. Our estimate of the power in 
the second harmonic, however, is larger than that in the first one, whereas they find 
the opposite. We attribute this difference to their use of a Gaussian distribution of 
sources in the Galactic disc rather than the (usually assumed) exponential that we 
adopted. Comparison between these two results suggests that it should be possible 
to infer the distribution of WD-WD binaries in our Galaxy by properly analyzing the 
harmonics of the variance of the galactic background measured by LISA. How this can 
be done will be the subject of a future work. 

As a next step in our analysis, we have estimated the cyclic spectra g r {f) fEa. ljMI 
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Cyclic spectra of WD-WD background 
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Figure 6. The main (k = 0) spectrum of simulated WD-WD background signal 
(red). 8 cyclic spectra (magenta) estimated from the simulated data. Spectral 
density of LISA instrumental noise (black) is shown for comparison. The 0th order 
spectrum contains LISA instrumental and hence it differs from the spectrum given 
in Figure HI 



with r = 1, 2, ...8) of the simulated WD-WD background signal. The number 8 comes 
from the theoretical predictions of the number of harmonics. The estimated absolute 
values of the cyclic spectra are shown in Figure where we also plot the spectrum 
of the LISA instrumental noise and the main spectrum (r = 0) estimated from the 
simulation. We find that the main spectrum and the first two cyclic spectra r = 1 and 
r = 2 have the largest magnitude and, over some frequency range, lie above the LISA 
instrumental noise. The remaining spectra are an order of magnitude smaller and are 
very noisy. We may also notice that all the cyclic spectra have roughly the same slope. 
This follows from the assumed independence between the location of the binaries in 
the Galaxy (D, A, (3) and their frequencies and chirp masses (M c ,u> s ). We also find 
the magnitude of the 2nd cyclic spectrum to be higher than the first, similarly to what 
we had for the harmonics of the variance. Note that we estimated the spectra from 
the time series consisting of the WD-WD background added to the LISA instrumental 
noise. 

Our analysis has shown that the LISA data will allow us to compute 17 
independent cyclic spectra (the 8 complex cyclic spectra g r (f),r = 1,2, ...8 and the 
real spectrum <7o(/)) of the WD-WD galactic background, 5 of which can be expected 
to be measured reliably. We have also shown that by performing generalized spectral 
analysis of the LISA data we will be able to derive more information about the WD- 
WD binary population (properties of the distribution of its parameters) than we would 
have by only looking at the ordinary go{f) spectrum. 



Simulation of the White Dwarf — White Dwarf galactic background in the LISA data.lA 



Acknowledgments 

The supercomputers used in this investigation were provided by funding from the Jet 
Propulsion Laboratory Institutional Computing and Information Services, and the 
NASA Directorates of Aeronautics Research, Science, Exploration Systems, and Space 
Operations. A.K. acknowledges support from the National Research Council under the 
Resident Research Associateship program at the Jet Propulsion Laboratory. This work 
was supported in part by Polish Science Committee Grant No. KBN 1 P03B 029 27. 
This research was performed at the Jet Propulsion Laboratory California Institute of 
Technology under contract with the National Aeronautics and Space Administration. 



References 



P. L. Bender, K. Danzmann, and the LISA Study Team, Laser Interferometer Space Antenna 
for the Detection of Gravitational Waves, Pre-Phase A Report, doc. MPQ 233 (Max-Planck- 
Institiit fur Quantenoptik, Garching, 1998). 

M.J. Benacquista, J. DeGoes, and D. Lunder, Class. Quantum Grav., 21, S509, (2004). 

D. Hils, P. L. Bender, and R. F. Webbink, Ap. J., 360, 75 (1900). 

Ch. Evans, Ickp Iben, and L. Smarr, Ap. J., 323, 129 (1987). 

N. Seto, Phys. Rev. D, 69, 123005 (2004). 

A. Krolak, M. Tinto, and M. Vallisneri, Phys. Rev. D, 70, 022003 (2004). 

J. A. Edlund, M. Tinto, A. Krolak, and G. Nelemans, In preparation. 

M. Tinto, F.B. Estabrook, and J.W. Armstrong, Phys. Rev. D 65, 082003 (2002). 

F. B. Estabrook, M. Tinto, and J.W. Armstrong, Phys. Rev. D, 62, 042002 (2000). 
D. Hils and P. L. Bender, Ap. J., 360, 75-94 (1990). 

D. Hils and P. L. Bender, Ap. J., 537, 334-341 (2000). 

P. L. Bender and D. Hils, Class. Quantum Grav., 14, 1439-1444 (1997). 

G. Nelemans, L. Ft. Yungelson, and S. F. Portegies-Zwart, A&A, 375, 890 (2001). 

G. Nelemans, L. R. Yungelson, and S. F. Portcgics-Zwart, Mon. Not. Roy. Astron. Soc, 349, 
181 (2004). 

Plots and data for Lisa sensitivity curve and WD-WD background due to Hils and Bender can 
be found on the web page ' http://www.srl. caltech.edu^shane/scnsitivity/MakeCurve. html ' . 
T. R. Marsh, G. Nelemans, and D. Steeghs, Mon. Not. Roy. Astron. Soc, 350, 113 (2004). 
J. W. Armstrong, F. B. Estabrook, and M. Tinto, Ap. J., 527, 814 (1999). 

G. Giampieri and A. G. Polnarev, Mon. Not. R. Astron. Soc, 291, 149-161 (1997). 

A.H. Nuttal, IEEE Transactions on Acoustics, Speech, and Signal Processing, ASSP-29, 84, 
(1981). 

H. L. Hurd, IEEE Transactions on Information Theory, 35, 350 (1989). 

D. B. Percival and A. T. Walden, Spectral Analysis for Physical Applications, Cambridge 
University Press (1993), Chapter 6. 



